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Abstract 

We discuss fitting hadronic Green functions versus time t to extract mass values in 
quenched lattice QCD. These data are themselves strongly correlated in t. With only a 
limited number of data samples, the method of minimising correlated is unreliable. We 
explore several methods of modelling the correlations among the data set by a few parame- 
ters which then give a stable and sensible fit even if the data sample is small. In particular 
these models give a reliable estimate of the goodness of fit. 



1 Introduction 

We assume that we have samples of unbiased estimators of quantities x{t) with t = 1 . . . D. 
Thus the data set is x^"'^ {t) where n = \ . . . N . We assume that the samples x^") (t) are statis- 
tically independent versus n for fixed t but may be correlated in t. Such a situation arises in 
lattice gauge theory calculations where there are independent configurations and D Green 
functions (vacuum expectation values of combinations of quark propagators and/or gauge links) 
are measured versus time separation t. An introduction to this topic in the context of lattice 
gauge theory is provided by Toussaint [|l|. 

The aim is to fit a given function F{t) which depends on P parameters Op. This function is 
to be fitted to the data samples x{t). Thus we require to find the following 

• the best values of the parameters ap. 

• the errors associated with these best fit parameters. 

• the confidence level that the fit represents the data sample. 
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A general discussion of this problem has been given |^] which has shown that with limited 
data samples N, it may be unrealistic to allow a general form for the D x D correlation matrix 
describing the correlations among the data at the D different t values. In particular the estimated 

was shown to be increased by a factor of 1 + (-D + 1) /N + 0{N~'^). This led to the conclusion 
that > 10(1) + 1) is needed for a reliable use of as a goodness of fit estimator in a correlated 

x' fit. 

Although this conclusion arose from studying a theoretical distribution which was uncorre- 
lated, the above increase is the same for any true distribution. Indeed the per degree 
of freedom will have the expected value (A^ — 1)/(A^ — D — 2). Furthermore the distribution 
of can be evaluated theoretically so that confidence levels can be derived. The major 
consequence of small sample size is that the eigenvalues of the correlation matrix are modified: 
in particular very small eigenvalues can arise. These correspond to a very narrow distribution 
in the corresponding eigen-direction and so can bias any fit considerably. 

In order to cope with this, it is necessary to use some extra theoretical input. In particular, 
one needs to assume a form of the correlation matrix or of its eigenvalue spectrum. This has 
been recognised before, and proposals have been made to truncate the eigenvalues according to 
the spirit of the SVD inverse of a singular matrix ^, ^j. We discuss this proposal and offer 
our suggestions for an improved treatment. We believe that our approach is more stable. 

Another way to attack this problem is to explore the theoretical expectations for the cor- 
relation matrix. In the case we are considering, the correlation between meson operators can 
be expressed in terms of vacuum expectation values of 4-quark operators [Q]. These can be 
estimated and this gives a motivation for a direct parametrization of the correlation matrix in 
terms of exponentials in |t — If the normalised correlation is given by one such exponential, 
then its inverse is a tridiagonal matrix. This is a very appealing model since it corresponds to a 
nearest neighbour linkage in the underlying probability distribution. This one-parameter model 
is very stable, but does not always give an accurate description of the data. A natural way to 
extend it is to consider a 5-diagonal inverse as an efficient parametrization. This 2-parameter 
model corresponds to a correlation matrix given by a specific combination of 2 exponentials and 
it fits the data well in many cases. 

We start by presenting typical data on the hadronic Green functions from quenched lattice 
studies. This enables us to analyse the behaviour of the correlation matrix which is at the heart 
of our exploration. We then compare various methods of modelling the correlation matrix and 
test their stability in correlated x^ fits. 
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2 Fits with correlated data 



The data x{t) we will consider are Green functions which are vacuum expectations of hadronic 
operators at times and t. Such data arc often referred to as hadron correlators but we will not 
use this description here since we wish to concentrate on a different correlation: that between 
the Green functions x{t) at different t values. 

The data sample themselves give a probability distribution 

1 ^ 

n=l 

We shall be interested in estimates of the probability distribution Q{X) of the averages X{t) of 
the data x{t). A smooth representation of Q{X) is needed for determining best fit parameters 
and to estimate the acceptability of such a fit. The natural parametrisation for Q is suggested 
by the central limit theorem. Provided that the underlying distributions of x{t) are sufficiently 
localised, then for large N, X{t) will be gaussianly distributed. We are specifically interested in 
the case where the different components x{t) are statistically correlated. Thus a general gaussian 
surface will be needed. 

Q{X) = HeM-^Ei^it) -m)M{t,t'){X{t') -X{t'))) 
t,t' 

with 

1 ^ 

x{t)=-j:x(-\t) 

n=l 

M{t,t') = NC-'^{t,t') ,where 
1 ^ _ _ 

n=l 

To find the best fit parameters then corresponds to maximising 

exp(-xV2) where x' = E(^(*' «) " '^{t))M{t, t'){F{t', a) - X{t')) 

t,t' 

with respect to for j3 = 1 . . . P. This is the usual correlated method. For sufficiently large 
this is a stable procedure and the expected value of is the number of degrees of freedom 
D — P. But for small N , the sampling fluctuation in C and hence M can give unreasonable fits. 
In order to avoid this bias, we aim to make a more stable model for C. 

Let us now study the properties of the correlation matrix C. It is a real symmetric positive- 
definite matrix for any number of samples N but it has rank A'^ — 1 so will have D — N +1 zero 
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eigenvalues if < D. In this latter case, its inverse is not defined. A commonly used prescription 
in such a case is the Singular Value Decomposition (SVD) inverse which corresponds to omitting 
the eigenmodes corresponding to the zero eigenvectors from the inverse. We will return to discuss 
the utility of this prescription. 

For many purposes, it is simpler to study the normalised correlation matrix which we define 



3 Correlations among hadronic operators 

In quenched lattice determinations of hadronic spectra and matrix elements, one studies vacuum 
expectation values of hadronic operators at times and t. Thus x{t) is the vacuum expectation 
value of hadronic operators <H{0)H(t) >. These determinations of x{t) are extracted from the 
quark propagators derived from inverting the fermion matrix in the given gauge field sample. 
This propagator inversion is very demanding computationally and is usually only evaluated 
for one source point (0,0). This implies that quark propagators to all values of t are equally 
sensitive to the region around this fixed source site. Thus between different gauge (vacuum) 
samples, all propagators will tend to be large/small as this fixed region is conducive/resistant to 
quark propagation. This argument shows that very prominent correlations are expected between 
hadron Green functions to different t values. We analyse some data to substantiate this. 

We present results for the normalised correlation C{t,t') in fig 1. Here the data come from 
a study of the vr and p mesons using local hadronic operators from a point source at (0, 0) to 
(y, t) summed over y to give a zero-momentum observable. There are = 60 independent 
configurations using a hopping parameter K = 0.14262 with the clover improved action |^]. The 
lattice has size 24^ x 48 at /3 = 6.2 and, for orientation, m^a ~ 0.17 at this hopping parameter 
value. We use t values 5 to 24 for this study. 

We find that C decreases with \t — t'\ but is relatively insensitive to t + t'. The decrease 
with |t — t'l is illustrated in fig 1. An exponential behaviour versus |i — i'| is expected from an 
analysis in terms of hadronic operators. 

Let us summarise this argument. Consider the case where x{t) is the vacuum expectation 
value of hadronic operators <H(0)H{t)>. Then 



The first term then will have contributions between t and t' from intermediate states of lowest 
energy m with the quantum numbers created by H, the same as those in x{t) itself. Between 
and t (for t' > t), the intermediate state will have the quantum numbers of H{0)H(0) and 



as 




C{t,t') =<H{0)H{0)H{t)H{t')> - <H{0)H{t)><H{0)H{t')> 
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so may have a lower energy which we write as 2M. Ignoring for the moment the disconnected 
term, and assuming that one state only dominates in each case, then 

For the case of the p for example: m = nip and M = m-j^ since a 27r state can couple to HH 
(ie to pp). Thus a small exponential rate of decrease is to be expected with exponent m,p — m.^. 
In this case the disconnected term in C will be relatively unimportant since it decreases faster 
than the connected term by exp(— 2(m — M)t). The contribution of excited intermediate states 
will modify this simple exponential behaviour except at large \t — t'\ where the lowest state 
dominates. The curve corresponding to the ground state exponential exp( — (rrip — m.j^)\t — is 
shown in fig.l, where it is seen to be a reasonable guide to the large \t — t'\ behaviour. 

For the case where a: is a pion observable, then both m and M arc and the correlation 
C{t,t') would be constant versus \t — t'\. In this case the disconnected term will be relatively 
important. The disconnected part will reduce the magnitude of C especially when relatively 
more disconnected terms are present - such as with a source summed over spatial position. 
When the source is at a fixed lattice position (0,0), as above, then the disconnected parts will 
cancel less completely and we expect C to remain large for large |i — Indeed it is larger for 
TT than p correlations as shown in fig. 1. 

For baryon spectra, the nucleon is the lowest lying 3 quark state and so C will be like the 
pion case above, while the A will behave analogously to the p above. 

4 Exponential models for correlations among the data 

Since there is some theoretical justification for an exponential decrease of C {t, t') with |t — 1'| , we 
consider first a simple and robust model with just one exponential. Since C(t,t') is normalised, 
this results in a one parameter model with parameter a 

(7(4,*') = e-"l*-*'l 

In practice the behaviour of C{t, t') is not exactly exponential, so one must choose a reference 
value of \t — t'\ = ta to determine a suitable value of a. We have found that using = 4 is a 
good choice. Then a is determined by averaging the D — ta values of C{t,t + ta) obtained from 
the sample data. This averaging along the off-diagonal of C also helps to reduce the sampling 
fluctuations. 

The test of the suitability of a model of C{t, t') for our purposes is that its inverse M is stable 
under fluctuations in the data sample used to model C{t,t'). Thus it is appropriate to study 
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the inverse of C{t,t'). For the case of a single exponential, the inverse is particularly simple: it 
is tridiagonal. It is given exactly in terms of a = exp(— a) by 

1 + ~ —a ~ 1 

M{t,t) = - 2' M{t,t±l) = - Mitm,tm) = - 2 

1 — 1 — 1 — 

where tm is the maximum or minimum t value in the matrix being inverted. 

Thus as a ^ 0, the elements of M{t,t') increase as a~^. In the limit of large D, one can 
estimate the smallest eigenvalue of C{t,t') (largest of M{t,t')) which is a/2 (2/a respectively) 
as a — > 0. Thus if a is reasonably well determined by the correlated data sample, then the value 
of M{t, t') will be stable under sample fluctuations. Thus the resultant fits will be stable too. 
This method provides a stable one parameter model for C{t,t'). One drawback of the model is 
that it does not reproduce very accurately any non-exponential behaviour of C{t,t'). This can 
be taken into account in a straightforward way by considering 2-exponential models for C{t,t'). 

Rather than consider an arbitrary 2-exponential model, we generalise the tri-diagonal feature 
of M(t, t') and look for 5-diagonal models instead. The algebra is now somewhat messier, but 
the conclusion is that a 5-diagonal model for M(i, t') corresponds to a particular 2-exponential 
model for C{t,t') with: 

^ (l-a2>il*-*'l+^-(l-«i>2l*-*'l+^ 
(ai — q;2)(1 + aia2) 
where to have a sensible interpretation we require < Reai, Rea2 < 1. 

If ai and a2 are complex, they must be complex conjugates. In this case the behaviour of 
C{t,t') versus \t — t'\ will be a damped oscillation. Although such a behaviour is not strictly 
excluded, it seems unreasonable to have anti-correlation at larger t values so we choose not to 
allow that possibility. Henceforward, we take cti and a2 to be real. 

Thus wc obtain ai and a2 by comparing the above expression for C{t,t') with the sample 
data. The two parameters can, for instance, be determined by making a least squares fit to 
C{t,t') at \t — t'\ = 1 and ta- A fit is needed because in some cases the data may not be 
reproducible exactly by the expression. We find that this 2-parameter assignment to the sample 
correlation is stable when inverted for use in fits. This is plausible since a 5-diagonal form of 
the inverse avoids very small eigenvalues of C{t,t'). 

We give the explicit formula for the inverse of the 5-diagonal matrix in the Appendix. 

5 Eigenvalue smoothing of the correlation matrix 

The essence of the problem is that sample values of C {t, t') may have very small eigenvalues 
and these influence unreasonably the inverse M(t, t') used in modelling the distribution of the 
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data. An obvious way to proceed is to modify these unreasonable eigenvalues by hand. One 
suggestion Q is to remove the largest eigenvalues of C{t, t') since they will have least influence on 
M{t, t'). This seems hard to justify and later suggestions |^ have been to remove the smallest 
eigenvalues of C{t,t'). This latter suggestion is in the spirit of the SVD inverse of a singular 
matrix: only the contributions from the non-zero eigenvalues are retained in the inverse. This 
eigenvalue truncation is clearly a rather brutal approximation to M{t, t'): its largest components 
are being removed. Indeed the gaussian surface modelling the probability distribution will be 
unconstrained in the direction of the deleted eigenmodes. A physical argument, for why this 
may be acceptable in practice, can be based on the observation |p that the smallest eigenvalues 
of C{t,t') usually correspond to eigenvectors which alternate in sign (versus t) and so are not 
very relevant to smooth fit functions. 

As we have argued, with a small sample size N, the D eigenvalues of the sample correlation 
matrix will be changed from their true values. The largest relative effect will come when there 
are several true eigenvalues of similar size - since those eigen directions mix fully. We find that 
the eigenvalue spectrum is densest at small eigenvalues. This again leads to the conclusion that 
the smallest eigenvalues of the sample correlation matrix are the most strongly affected. 

The SVD approach replaces the smallest (or zero) eigenvalues by very large values (since this 
corresponds to ignoring those terms in the inverse). We have explored this situation and found 
a less discontinuous way to modify the smallest eigenvalues Aj of C{t,t'). We keep the largest 
E eigenvalues substantially unchanged and then rearrange the remaining smaller eigenvalues to 
avoid extremely small ones. Following from the motivation that like eigenvalues mix most, we 
propose to redistribute the small eigenvalues by replacing them by their average to avoid very 
small values. 

After some experimentation, we propose the following explicit scheme to replace the D 
eigenvalues Aj of the sample normalised correlation C{t,t') (with convention Aj > Xi+i) with 
new eigenvalues A^: 

X'i = K Max(Aj, Xmin) 

I D _ 1 ^ 

where Xmin = t\ ^ Yl ^ ^ = T^Yl Max(Ai,A™„) 

The eigenvectors of C{t,t') and thus of its inverse M{t,t') are retained unchanged. Thus our 
procedure removes any very small eigenvalues of C{t,t') and replaces them with the average of 
the D — E smallest eigenvalues while retaining the property that the trace of C{t,t') is D. The 
procedure also ensures a smooth eigenvalue distribution by allowing eigenvalues larger than this 
average to be retained. 

We tested this assignment with the same data as used above. Of course foic E = D — 1, 
the method is equivalent to an exact inversion of the sample correlation matrix. For E > 2, 
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this method provides a stable model of the correlation matrix from the sample data. As E is 
increased, the model reproduces more exactly the sample correlation matrix - but at the expense 
of including unreasonable fluctuations if the sample size is too small. A compromise is to retain 
E ~ \/iV exact eigenvalues when there are A'^ samples. 

6 Practical test of models 

Here we apply the various models described above to some typical real data. Since we need a 
large number of samples to give an accurate data set we chose some APE data on p meson 
Green functions x{t) = < H(0)H(t) >. Here local-source and local-sink operators H are used 
for creating and destroying a p meson. The data set has A^ = 420 samples of x{t) with the clover 
fermionic action at K = .14190 on a 18^ X 64 lattice at /? = 6.2. 

The data sample is large enough to allow a full correlated fit to the observed x{t). An 
acceptable fit (x^/d.o.f. = 11.11/9) to x{t) with one exponential (actually a cosh is used) is 
found for the t-range 14 to 24. This 2-parameter fit corresponds to requiring a plateau in the 
effective mass. The eigenvalues of the normalised correlation matrix for this data set are shown 
in fig 2 by the continuous line. 

Here our intention is not to obtain the most precise values of nip but to illustrate the stability 
of various fitting prescriptions. Thus we take this full data sample as the true result and explore 
fits using smaller subsets of the 420 data. For example we take 14 blocks of 30 data. For each 
such set of 30, we perform a fit to the same t-range. The intention is to check whether the 
X^/d.o.f. is similar to the true value from the full data set. As well as the average value of 
X^/d.o.f., one may study its distribution so that confidence levels can be extracted. We do not 
pursue this here. 

A valid criterion for goodness of fit is of importance. The usual method is to use to decide 
if a fit is acceptable over a given t-range. Thus an erroneous value will lead to an increase or 
decrease in the t-range chosen as acceptable. This in turn will bias the fitted parameters such 
as nip. For example, an uncorrelated fit will yield much lower which will then suggest lower 
t values being included. This will tend to increase nip since the effective mass is a decreasing 
function of t in this case. 

Returning now to the comparison: fig 3 shows the results of the average of 14 fits to different 
samples of size A^ = 30. The full correlated fit (n = 11) has a higher on average by 40% 
than the 'true' result. This is entirely expected from the factor 1 + (Z? + 1)/A^ = 1.39 increase 
in x"^ predicted |p. The eigenvalues from one such sample fit to A^ = 30 are shown in fig 2. 
Here the phenomenon of very small eigenvalues appearing for small sample size is clear. Thus 
a direct use of a correlated fit to A^ = 30 samples is likely to be biassed by those spurious 
small eigenvalues. We consider now some of the models introduced above to counter this while 
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retaining a reasonable estimate of the goodness of fit. 

As shown in fig 3, the uncorrelated fit has a very much reduced ^ as expected. The 1- and 
2-exponential models of section 4 do much better in estimating x^- The eigenvalue smoothing 
model of section 5 also does well when 4 to 8 eigenvalues are retained exactly rather than 
all eigenvalues (n = 11). Note that this is indeed consistent with our previous estimate that 
approximately \/iV exact eigenvalues can be relied on. The modifications to the eigenvalues 
from the smoothing model are shown in fig 2 for one sample and 6 exact eigenvalues. Indeed 
the modification does alter the eigenvalues from the sample in the direction of the true values. 

Also shown in fig 3 are the results using an SVD definition of the inverse of C(t, t') in which 
n eigenvalues are retained exactly and the remainder are discarded. It is possible to estimate the 
expected value of in this approach yielding N{n — 2)/{N — n — 2) which is quite close to the 
values in fig 3. Thus a corrected estimator of the goodnes of fit is essential in this approach. The 
eigenvalue distribution is not smooth since the deletion of eigenmodes is equivalent to replacing 
the D — n smallest eigenvalues in fig 2 by infinite values. 

For the samples of = 30, the situation is that the correlations appear stronger than the true 
distributions. This can be understood from the earlier discussion of the eigenvalues of C{t,t'). 
The true correlation matrix C has 11 eigenvalues with values in the range 0.00438 to 9.82 whereas 
from subsets of 30 samples the smallest eigenvalue fiuctuates between 0.00052 to 0.002844. These 
reductions in the magnitude of the smallest eigenvalues correspond to narrower probability 
distributions and hence stronger correlations. The implications of this for the goodness of fit 
have been discussed, but we also need to check that the fitting procedures do not upset the 
fitted parameters themselves. For our chosen i-range the correlated fit to the full 420 samples 
yields rupa = 0.351(5). A range of other fits to the same full sample (ie uncorrelated, 5-diagonal, 
smoothed, etc) give essentially the same value of rrip. Thus it is only the goodness of fit that 
depends on the correlation model used. For the fits to subsets of 30 samples, we find that an 
uncorrelated fit does give the same result when averaged over the 14 independent blocks. The 
various correlation models all give a result somewhat higher (nipa ~ 0.358). Thus we have some 
evidence that the uncorrelated fit is the most stable for determining the fit parameters when 
the sample size is reduced. 



7 Conclusion 



Data appropriate to hadron propagation in lattice gauge theory calculations are very strongly 
correlated. A full correlated fit to such data can be biassed unless the sample size (N) is suffi- 
ciently large compared with the number of data points (D). The main effect is the appearance 
of spurious small eigenvalues of the correlation matrix. These increase the correlation in the 
sample. This increases by a factor of {N — 1)/ {N — D — 2). It can also bias the fit parameters. 
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Here we propose models which can amchorate this. Such models depend, to some extent, on 
an understanding of the expected form of the correlation matrix. Thus such methods arc not 
completely general. For the applications considered here, we find two promising avenues. One is 
to require that the normalised correlation matrix has a tri-diagonal or 5-diagonal inverse. The 
other is to average the smallest eigenvalues so that spurious small values are removed. Both of 
these models give reasonable estimates of the goodness of fit even with quite small sample size. 
The goodness of fit is important because it determines the range of data (eg. the i-range) to be 
fitted. This then will influence the fitted parameters in turn. 

For the actual best determination of the fitted parameters for a small sample, we find that 
an uncorrelated fit is stable and thus an attractive proposition. 

In fitting to hadron Green functions, it is preferable to make a simultaneous fit to several 
operators. For example smeared or fuzzed operators can be used at either sink or source as well 
as local ones. The observables for different operators will be correlated in general. We have 
discussed the correlation in t, but this correlation among different operators will need somewhat 
different models. The most attractive route is to use the method of smoothing the eigenvalues 
of the correlation matrix which is now treated as a PT x PT matrix if there are T t-values and 
P operators. 
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A Appendix 

Consider the matrix given by 

_ (l-a2^)«il^-^'l+^-(l-ai^)«2l^-^'l+^ 

0(t,t j — r— 

(«! - Q!2)(l + aia2) 

Its inverse M{t, t') is 5-diagonal and is given by the following expressions in terms of p\ 
C{t, t ± 1) and p2 = C{t, t ± 2) with d=l- 2pi^ + 2pi^p2 - 



M{t,t±l) = — ^ , M{t,t±2) = P^—Pl , M{t,t) = l-2piM{t,t±l)-2p2M{t,t±2) , 



dil-pi) 



d 
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M{tm, tm) = , M(t„, tn) = " PlM(t, t ± 1) , Af(t„, t„ ± 1) = "^^^^J^- 

where tm is the maximum or minimum value of t and tn is the next to maximum or minimum. 
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Figure 1: The normalised correlation C{t,t') versus the relative t difference for pion and rho 
Green functions from the data of section 3. The continuous line is the 2-exponential model, 
the dotted line is one exponential and the dot-dashed curve is the expected behaviour from 
theoretical analysis of the 4 point function. 
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Figure 2: The eigenvalues of the 11 x 11 normahsed correlation C{t, t') for the data described in 
section 6. The full line is the result from N = 420 configurations. The sample values (diamonds) 
are from a subset of A?^ = 30. The squares represent the smoothed prescription of section 5 with 
6 exact eigenvalues. 
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Figure 3: The for the P = 2 parameter correlated fit to -D = 11 data. The fuU hne is the 
result for the correlated fit to the full data set of N = 420. The expected is given by the 
number of degrees of freedom D — P = 9. The same fit was made to samples of = 30 data. 
The average and standard deviation of 14 such independent fits are shown. The uncorrclatcd, 
1-exponential, and 2-exponential models of section 4 are shown by diamonds. The eigenvalue 
smoothed models of section 5 are shown by squares. For comparison an SVD-truncated model 
is also shown (bursts). 
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